function [Connectedness] = ComputeSimilarity_Polk(X)

%%%% This function computes the Anton-Polk (2014) similarity measure.

X = sparse(X);
B = sparse(double(X>0));
TNA = sum(X')';

% reduce number of columns
i = find(sum(B)<=1);
X(:,i) =[];
B(:,i) =[];

[N,K] = size(X);

%%%%%%%% compute connectedness

% Connectedness = (X*B' + B*X');

% split up the matrix multiplication to speed up the process
Connectedness = zeros(N);
nsteps = 5;
nn = floor(K/nsteps);

for i = 1:nsteps
    if i < nsteps
    tmp_X = X(:,(i-1)*nn+1:i*nn);
    tmp_B = B(:,(i-1)*nn+1:i*nn);
    else
    tmp_X = X(:,(i-1)*nn+1:end);
    tmp_B = B(:,(i-1)*nn+1:end);    
    end
    Connectedness = Connectedness + (tmp_X*tmp_B' + tmp_B*tmp_X');
end

Connectedness = Connectedness./(repmat(TNA,1,N) + repmat(TNA,1,N)');

for i = 1:N
    Connectedness(i,i) = NaN;
end

end